home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C/C++ Users Group Library 1996 July
/
C-C++ Users Group Library July 1996.iso
/
listings
/
v_12_06
/
gage
/
divdif.c
< prev
Wrap
C/C++ Source or Header
|
1994-02-06
|
5KB
|
191 lines
/* Divided Difference Interpolation */
/* Copyright 1994 by Philip Gage */
#include <stdio.h>
#include <math.h>
#define POINTS 10 /* Max number of data points */
/* Function prototypes */
void interpolate(int,int,int,double*,double*,double*);
double evaluate (int, int, double*, double);
void print_function (int, int, double*);
void print_table(int,int,int,double*,double*,double*);
/* Generate polynomial or exponential function
polynomial 1=polynomial, 0=exponential
degree Degree of desired solution
n Number of data points
x,y Data point values
coefficient Output coefficient array */
void interpolate (int polynomial, int degree, int n,
double *x, double *y, double *coefficients)
{
double xpower, data[POINTS], diff[POINTS];
int i, j, k;
/* Copy input y values to data array */
for (i = 0; i < n; i++)
data[i] = y[i];
/* Initialize coefficients to identity element */
for (i = 0; i <= degree; i++)
if (polynomial)
coefficients[i] = 0.0;
else
coefficients[i] = 1.0;
/* Solve each high order term */
for (i = degree; i >= 0; i--) {
/* Copy data values to difference array */
for (j = 0; j < n; j++)
diff[j] = data[j];
/* Compute divided differences or quotients */
for (j = 1; j <= i; j++)
for (k = 0; k < n-j; k++)
if (polynomial)
diff[k] = (diff[k+1] - diff[k]) /
(x[k+j] - x[k]);
else
diff[k] = pow(diff[k+1] / diff[k],
1.0/(x[k+j] - x[k]));
/* Average current difference row */
coefficients[i] = 0.0;
for (j = 0; j < n-i; j++)
coefficients[i] += diff[j];
coefficients[i] /= n-i;
/* Remove high order term from data */
for (j = 0; j < n; j++) {
/* Compute x[j] raised to i power */
xpower = 1.0;
for (k = 0; k < i; k++)
xpower *= x[j];
if (polynomial)
data[j] -= coefficients[i] * xpower;
else
data[j] /= pow(coefficients[i],xpower);
}
}
}
/* Evaluate function at xvalue using Horner's rule */
double evaluate (int polynomial, int degree,
double *coefficients, double xvalue)
{
double yvalue;
int i;
if (polynomial)
yvalue = 0.0;
else
yvalue = 1.0;
for (i = degree; i >= 0; i--)
if (polynomial)
yvalue = yvalue * xvalue + coefficients[i];
else
yvalue = pow(yvalue,xvalue) * coefficients[i];
return yvalue;
}
/* Print symbolic function terms using ^ for power */
void print_function (int polynomial, int degree,
double *coefficients)
{
int i, identity, printing = 0;
char operator1, operator2;
/* Setup for polynomial or exponential */
if (polynomial) {
identity = 0.0;
operator1 = '+';
operator2 = '*';
}
else {
identity = 1.0;
operator1 = '*';
operator2 = '^';
}
printf("\nF(X) = ");
for (i = degree; i >= 0; i--) {
if (coefficients[i] != identity) {
if (printing)
printf(" %c ",operator1);
printf("%g",coefficients[i]);
if (i > 0)
printf("%cX",operator2);
if (i > 1)
printf("^%d",i);
printing = 1;
}
}
printf("\n");
}
/* Print table of results */
void print_table (int polynomial, int degree, int n,
double *x, double *y, double *coefficients)
{
int i;
double result, error, percent;
printf("%15s%15s%15s%15s%15s %\n",
"X","Y","F(X)","Error","Error");
for (i = 0; i < n; i++) {
result = evaluate(polynomial,degree,
coefficients,x[i]);
error = result - y[i];
if (y[i] != 0.0)
percent = 100.0 * error / y[i];
printf("%15g%15g%15g%15g%15g %\n\n",
x[i],y[i],result,error,percent);
}
}
void main()
{
int i, n, degree, polynomial, choice;
double x[POINTS], y[POINTS], coefficients[POINTS];
/* Read number of data points from 2 to POINTS */
do {
printf("Number of data points: ");
scanf("%d",&n);
} while (n<2 || n>POINTS);
/* Read each data point from user */
for (i = 0; i < n; i++) {
printf("x%d: ",i);
scanf("%lf",&x[i]);
printf("y%d: ",i);
scanf("%lf",&y[i]);
}
/* Generate solutions until Quit selected */
while (1) {
printf("1 Polynomial interpolation\n");
printf("2 Exponential interpolation\n");
printf("3 Quit program\n");
printf("Choice: ");
scanf("%d",&choice);
if (choice < 1 || choice > 2) break;
polynomial = (choice == 1);
/* Read degree of solution from 1 to n-1 */
do {
printf("Degree: ");
scanf("%d",°ree);
} while (degree < 0 || degree >= n);
/* Compute function and print results */
interpolate(polynomial,degree,n,x,y,coefficients);
print_function(polynomial, degree, coefficients);
print_table(polynomial,degree,n,x,y,coefficients);
}
}